##################################################################################################################
# Code for Results in Appendix Section A4
##################################################################################################################
write.table(data.frame(Description=c("","","","Appendix Section A4: Number and Percentage of Excluded Observations")),row.names=F,col.names=F,quote=F)

##################################################################################################################
# To run, please set the working directory to match the location of the folder containing all replication files.
##################################################################################################################
#setwd("Replication_PSRM_ST24")

#if not yet created, create "output" folder
if("output"%in%list.files()==F){dir.create("output")}

#Clear environment
rm(list=ls());gc();gc();gc();gc()

#If not yet installed, install package
if(length(find.package("dplyr",quiet=T))==0){install.packages("dplyr")}

###packages
library(dplyr)
#packageVersion("dplyr") #[1] '1.1.4'

#Set POSIX locale
Sys.setlocale(locale="C")

##################
#Load main dataset
##################
loball5ciA<-readRDS("Data/PSRM_main_dataset.rds")


####################
# Prepare Analyes
####################

#####
#within dyad deviations

#contract and in-house
agg67c<-loball5ciA %>%
  filter(type%in%c("contract","inhouse")==T) %>%
  group_by(principal_cleaned,lobname_cleaned, halfyear,type,session,both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>% 
  summarise(.groups="keep",sum_hours=sum(hours,na.rm=T),sum_amount=sum(Rcomparable_total_amounts_cleaned_exp,na.rm=T))

agg67d2<-agg67c
#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg67d2$sum_hours_cr0<-agg67d2$sum_hours
agg67d2$sum_hours_cr0[agg67d2$sum_hours_cr0>0]<-agg67d2$sum_hours_cr0[agg67d2$sum_hours_cr0>0]/min(agg67d2$sum_hours_cr0[agg67d2$sum_hours_cr0>0])
agg67d2$sum_hours_cr0[agg67d2$sum_hours_cr0==0]<-1

length(which(agg67d2$sum_hours==0&agg67d2$sum_amount==0)) #9207

write.table(data.frame(Description=c("","Results:")),row.names=F,col.names=F,quote=F)
write.table(data.frame(Description=c("","Footnote 16 (and Appendix A4): Percentage of observations with zero expenditures that have non-zero hours associated with them:")),row.names=F,col.names=F,quote=F)
print(round(length(which(agg67d2$sum_hours==0&agg67d2$sum_amount==0))/length(which(agg67d2$sum_amount==0)),2)) #90%

#count number of observations
n671<-nrow(agg67d2) #[1] 43051
table(agg67d2$type) 
# contract  inhouse 
# 27775    15276

#create variables for analyses
agg67d2 <-agg67d2 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned,type,session) %>%
  mutate(
         dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
         dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
         )

n672<-nrow(agg67d2) #[1] 32813
#43051-32813=10238
#10238/43051=0.24
table(agg67d2$type) 
# contract  in-house 
# 20254    12559

agg67d2$count<-1
agg67d2 <-agg67d2 %>% 
  group_by(principal_cleaned,lobname_cleaned,type,session) %>%
  mutate(num_dyad_obs=sum(count)) 
#exclude observations when num_dyad_obs==1
agg67d2<-agg67d2[which(agg67d2$num_dyad_obs>1),]
#count obs
n673<-nrow(agg67d2) #[1] 30888
#32813-30888=1925
#1925/43051=0.04
table(agg67d2$type) 
# contract  inhouse 
# 18890    11998

#12163/43051=28% total removed (10238+1925)

###contract, half-yearly
agg6c<-loball5ciA %>%
  filter(type%in%c("contract")==T) %>%
  group_by(principal_cleaned,lobname_cleaned, halfyear,session,both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>% 
  summarise(.groups="keep",sum_hours=sum(hours,na.rm=T),sum_amount=sum(Rcomparable_total_amounts_cleaned_exp,na.rm=T))

#get total amount of compensation (by dyad) in biennium
agg6c <- agg6c %>% 
  group_by(lobname_cleaned,halfyear) %>%
  mutate(count_client=length(unique(principal_cleaned)))

agg6d2<-agg6c
#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg6d2$sum_hours_cr0<-agg6d2$sum_hours
agg6d2$sum_hours_cr0[agg6d2$sum_hours_cr0>0]<-agg6d2$sum_hours_cr0[agg6d2$sum_hours_cr0>0]/min(agg6d2$sum_hours_cr0[agg6d2$sum_hours_cr0>0])
agg6d2$sum_hours_cr0[agg6d2$sum_hours_cr0==0]<-1

n61<-nrow(agg6d2) # 27775

agg6d2 <-agg6d2 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned,session) %>%
  mutate(
         dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
         dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
         )

#count #obs
n62<-nrow(agg6d2) #20254
#27775-20254=7521
#7521/27775=0.27
agg6d2$count<-1
agg6d2 <-agg6d2 %>% 
  group_by(principal_cleaned,lobname_cleaned,session) %>%
  mutate(num_dyad_obs=sum(count)) 
#exclude observations when num_dyad_obs==1
agg6d2<-agg6d2[which(agg6d2$num_dyad_obs>1),]
#count obs
n63<-nrow(agg6d2) #18890
#20254-18890= 1364
# 1364/27775=0.05

#####in-house

#in-house
agg7c<-loball5ciA %>%
  filter(type%in%c("inhouse")==T)%>%
  group_by(principal_cleaned,lobname_cleaned,halfyear,session,both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>% 
  summarise(.groups="keep",sum_hours=sum(hours,na.rm=T),sum_amount=sum(Rcomparable_total_amounts_cleaned_exp,na.rm=T))

#get number of clients for whom lobbyist is working in halfyear
agg7c <- agg7c %>% group_by(lobname_cleaned,halfyear) %>%
  mutate(count_client=length(unique(principal_cleaned)))



agg7d2<-agg7c
#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg7d2$sum_hours_cr0<-agg7d2$sum_hours
agg7d2$sum_hours_cr0[agg7d2$sum_hours_cr0>0]<-agg7d2$sum_hours_cr0[agg7d2$sum_hours_cr0>0]/min(agg7d2$sum_hours_cr0[agg7d2$sum_hours_cr0>0])
agg7d2$sum_hours_cr0[agg7d2$sum_hours_cr0==0]<-1

n71<-nrow(agg7d2) #15276

#create variables for analyses
agg7d2 <-agg7d2 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned,session) %>%
  mutate(
         dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
         dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
         )

n72<-nrow(agg7d2) # 12559
# 15276-12559= 2717
# 2717/15276=0.18

agg7d2$count<-1
agg7d2 <-agg7d2 %>% group_by(principal_cleaned,lobname_cleaned,session) %>%
  mutate(num_dyad_obs=sum(count)) 
#exclude observations when num_dyad_obs==1
agg7d2<-agg7d2[which(agg7d2$num_dyad_obs>1),]

n73<-nrow(agg7d2) # 11998
# 12559-11998=561
# 561/15276=0.04


############################################################
############ Aggregated to yearly
############################################################

agg67d3<-agg67c %>%
  group_by(principal_cleaned,lobname_cleaned,year=substr(halfyear,1,4),session, type,both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>%
  summarise(.groups="keep",sum_hours=sum(sum_hours,na.rm=T),sum_amount=sum(sum_amount,na.rm=T))

#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg67d3$sum_hours_cr0<-agg67d3$sum_hours
agg67d3$sum_hours_cr0[agg67d3$sum_hours_cr0>0]<-agg67d3$sum_hours_cr0[agg67d3$sum_hours_cr0>0]/min(agg67d3$sum_hours_cr0[agg67d3$sum_hours_cr0>0])
agg67d3$sum_hours_cr0[agg67d3$sum_hours_cr0==0]<-1

n67a1<-nrow(agg67d3) #22611

#create variables for analyses
agg67d3 <-agg67d3 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(
         dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
         dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
         )

n67a2<-nrow(agg67d3) #18902
#22611-18902=3709
#3709/22611=0.16
agg67d3$count<-1

agg67d3 <-agg67d3 %>% 
  group_by(principal_cleaned,lobname_cleaned,type) %>%
  mutate(num_dyad_obs=sum(count)) 
#exclude observations when num_dyad_obs==1
agg67d3<-agg67d3[which(agg67d3$num_dyad_obs>1),]

n67a3<-nrow(agg67d3) #17404
#18902-17404=1498
#1498/22611=0.07

agg67d3<-agg67d3[order(agg67d3$principal_cleaned,agg67d3$lobname_cleaned,as.numeric(agg67d3$year),agg67d3$type),]
agg67d3<-agg67d3 %>% group_by(principal_cleaned,lobname_cleaned,type) %>%
  mutate(nyears=cumsum(count))

###contract lobbyists
agg6d3<-agg6c %>%
  group_by(principal_cleaned,lobname_cleaned,year=substr(halfyear,1,4),session, both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>%
  summarise(.groups="keep",sum_hours=sum(sum_hours,na.rm=T),sum_amount=sum(sum_amount,na.rm=T))

n6a1<-nrow(agg6d3) #14594

#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg6d3$sum_hours_cr0<-agg6d3$sum_hours
agg6d3$sum_hours_cr0[agg6d3$sum_hours_cr0>0]<-agg6d3$sum_hours_cr0[agg6d3$sum_hours_cr0>0]/min(agg6d3$sum_hours_cr0[agg6d3$sum_hours_cr0>0])
agg6d3$sum_hours_cr0[agg6d3$sum_hours_cr0==0]<-1

#create variables for analyses
agg6d3 <-agg6d3 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(
    dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
    dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
    )

n6a2<-nrow(agg6d3) # 11876
#14594-11876=2718
#2718/14594=0.11
agg6d3$count<-1

agg6d3 <-agg6d3 %>% group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(num_dyad_obs=sum(count)) 

#exclude observations when num_dyad_obs==1
agg6d3<-agg6d3[which(agg6d3$num_dyad_obs>1),]
n6a3<-nrow(agg6d3) # 10751
#11876-10751=1125
#1125/14594=0.08

agg6d3<-agg6d3[order(agg6d3$principal_cleaned,agg6d3$lobname_cleaned,as.numeric(agg6d3$year)),]
agg6d3<-agg6d3 %>% group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(nyears=cumsum(count))

#
agg7d3<-agg7c %>%
  group_by(principal_cleaned,lobname_cleaned,year=substr(halfyear,1,4),session, both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>%
  summarise(.groups="keep", sum_hours=sum(sum_hours,na.rm=T),sum_amount=sum(sum_amount,na.rm=T))

n7a1<-nrow(agg7d3) # 8017

#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg7d3$sum_hours_cr0<-agg7d3$sum_hours
agg7d3$sum_hours_cr0[agg7d3$sum_hours_cr0>0]<-agg7d3$sum_hours_cr0[agg7d3$sum_hours_cr0>0]/min(agg7d3$sum_hours_cr0[agg7d3$sum_hours_cr0>0])
agg7d3$sum_hours_cr0[agg7d3$sum_hours_cr0==0]<-1

#create variables for analyses
agg7d3 <-agg7d3 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(
    dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
    dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
    )

n7a2<-nrow(agg7d3) #7026
#8017-7026=991
#991/8017=0.06
agg7d3$count<-1

agg7d3 <-agg7d3 %>% group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(num_dyad_obs=sum(count)) 

#exclude observations when num_dyad_obs==1
agg7d3<-agg7d3[which(agg7d3$num_dyad_obs>1),]
n7a3<-nrow(agg7d3) #6653
#7026-6653=373
#373/8017 5%

agg7d3<-agg7d3[order(agg7d3$principal_cleaned,agg7d3$lobname_cleaned,as.numeric(agg7d3$year)),]
agg7d3<-agg7d3 %>% group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(nyears=cumsum(count))

# png("output/PSRM_fig1.png",width=12,height=8.83,units="in", res=300)
#   par(mfrow=c(2,3))
#   #Panel 1
#   tt<-cor.test(agg67d2$dif_log_mean_amount,agg67d2$dif_log_mean_hours_cr0)
#   plot(main=c("(1) All Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5, x=agg67d2$dif_log_mean_amount,y=agg67d2$dif_log_mean_hours_cr0,
#        xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
#   text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
#   legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
#   lines(lowess(x=agg67d2$dif_log_mean_amount,y=agg67d2$dif_log_mean_hours_cr0,f=0.68),lty=2)
#   #Panel 2
#   tt<-cor.test(agg6d2$dif_log_mean_amount,agg6d2$dif_log_mean_hours_cr0)
#   plot(main=c("(2) Contract Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg6d2$dif_log_mean_amount,y=agg6d2$dif_log_mean_hours_cr0,
#        xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
#   text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
#   legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
#   lines(lowess(x=agg6d2$dif_log_mean_amount,y=agg6d2$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
#   #Panel 3
#   tt<-cor.test(agg7d2$dif_log_mean_amount,agg7d2$dif_log_mean_hours_cr0)
#   plot(main=c("(3) In-House Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg7d2$dif_log_mean_amount,y=agg7d2$dif_log_mean_hours_cr0,
#        xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
#   text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
#   legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
#   lines(lowess(x=agg7d2$dif_log_mean_amount,y=agg7d2$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
#   #Panel 4
#   tt<-cor.test(agg67d3$dif_log_mean_amount,agg67d3$dif_log_mean_hours_cr0)
#   plot(main=c("(4) All Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg67d3$dif_log_mean_amount,y=agg67d3$dif_log_mean_hours_cr0,
#        xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
#   text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
#   legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
#   lines(lowess(x=agg67d3$dif_log_mean_amount,y=agg67d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
#   #Panel 5
#   tt<-cor.test(agg6d3$dif_log_mean_amount,agg6d3$dif_log_mean_hours_cr0)
#   plot(main=c("(5) Contract Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg6d3$dif_log_mean_amount,y=agg6d3$dif_log_mean_hours_cr0,
#        xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
#   text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
#   legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
#   lines(lowess(x=agg6d3$dif_log_mean_amount,y=agg6d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
#   #Panel 6
#   tt<-cor.test(agg7d3$dif_log_mean_amount,agg7d3$dif_log_mean_hours_cr0)
#   plot(main=c("(6) In-House Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg7d3$dif_log_mean_amount,y=agg7d3$dif_log_mean_hours_cr0,
#        xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
#   text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
#   legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
#   lines(lowess(x=agg7d3$dif_log_mean_amount,y=agg7d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
# 
# dev.off()
# 
# write.table(data.frame(Description=c("(N.B. The file for this figure is saved as 'PSRM_fig1.png' in the 'output' folder.)")),row.names=F,col.names=F,quote=F)
# 
# #Include code again so that Figure A6 also shown in R if code file is run using source()
# par(mfrow=c(2,3))
# #Panel 1
# tt<-cor.test(agg67d2$dif_log_mean_amount,agg67d2$dif_log_mean_hours_cr0)
# plot(main=c("(1) All Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5, x=agg67d2$dif_log_mean_amount,y=agg67d2$dif_log_mean_hours_cr0,
#      xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
# text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
# legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
# lines(lowess(x=agg67d2$dif_log_mean_amount,y=agg67d2$dif_log_mean_hours_cr0,f=0.68),lty=2)
# #Panel 2
# tt<-cor.test(agg6d2$dif_log_mean_amount,agg6d2$dif_log_mean_hours_cr0)
# plot(main=c("(2) Contract Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg6d2$dif_log_mean_amount,y=agg6d2$dif_log_mean_hours_cr0,
#      xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
# text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
# legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
# lines(lowess(x=agg6d2$dif_log_mean_amount,y=agg6d2$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
# #Panel 3
# tt<-cor.test(agg7d2$dif_log_mean_amount,agg7d2$dif_log_mean_hours_cr0)
# plot(main=c("(3) In-House Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg7d2$dif_log_mean_amount,y=agg7d2$dif_log_mean_hours_cr0,
#      xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
# text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
# legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
# lines(lowess(x=agg7d2$dif_log_mean_amount,y=agg7d2$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
# #Panel 4
# tt<-cor.test(agg67d3$dif_log_mean_amount,agg67d3$dif_log_mean_hours_cr0)
# plot(main=c("(4) All Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg67d3$dif_log_mean_amount,y=agg67d3$dif_log_mean_hours_cr0,
#      xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
# text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
# legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
# lines(lowess(x=agg67d3$dif_log_mean_amount,y=agg67d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
# #Panel 5
# tt<-cor.test(agg6d3$dif_log_mean_amount,agg6d3$dif_log_mean_hours_cr0)
# plot(main=c("(5) Contract Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg6d3$dif_log_mean_amount,y=agg6d3$dif_log_mean_hours_cr0,
#      xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
# text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
# legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
# lines(lowess(x=agg6d3$dif_log_mean_amount,y=agg6d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
# #Panel 6
# tt<-cor.test(agg7d3$dif_log_mean_amount,agg7d3$dif_log_mean_hours_cr0)
# plot(main=c("(6) In-House Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg7d3$dif_log_mean_amount,y=agg7d3$dif_log_mean_hours_cr0,
#      xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
# text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
# legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
# lines(lowess(x=agg7d3$dif_log_mean_amount,y=agg7d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)



write.table(data.frame(Description=c("","Overall number of half-yearly observations before filtering:")),row.names=F,col.names=F,quote=F)
print(n671)
write.table(data.frame(Description=c("","Number of half-yearly observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(n671-n672)
write.table(data.frame(Description=c("","Percentage of half-yearly observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n671-n672)/n671,2))
write.table(data.frame(Description=c("","Number of half-yearly observations excluded due to dyads with only one half-year with non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(c(n672-n673))
write.table(data.frame(Description=c("","Percentage of half-yearly observations excluded due to dyads with only one half-year with non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n672-n673)/n671,2))
write.table(data.frame(Description=c("","Total percentage of half-yearly observations excluded:")),row.names=F,col.names=F,quote=F)
print(round(c(n671-n673)/n671,2))

write.table(data.frame(Description=c("","Number of half-yearly contract lobbyist observations before filtering:")),row.names=F,col.names=F,quote=F)
print(n61)
write.table(data.frame(Description=c("","Number of half-yearly contract lobbyist observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(n61-n62)
write.table(data.frame(Description=c("","Percentage of half-yearly contract lobbyist observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n61-n62)/n61,2))
write.table(data.frame(Description=c("","Number of half-yearly contract lobbyist observations excluded due to dyads with only one half-year with non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(c(n62-n63))
write.table(data.frame(Description=c("","Percentage of half-yearly contract lobbyist observations excluded due to dyads with only one half-year with non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n62-n63)/n61,2))
write.table(data.frame(Description=c("","Total percentage of half-yearly contract lobbyist observations excluded:")),row.names=F,col.names=F,quote=F)
print(round(c(n61-n63)/n61,2))

write.table(data.frame(Description=c("","Number of half-yearly in-house lobbyist observations before filtering:")),row.names=F,col.names=F,quote=F)
print(n71)
write.table(data.frame(Description=c("","Number of half-yearly in-house lobbyist observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(n71-n72)
write.table(data.frame(Description=c("","Percentage of half-yearly in-house lobbyist observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n71-n72)/n71,2))
write.table(data.frame(Description=c("","Number of half-yearly in-house lobbyist observations excluded due to dyads with only one half-year with non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(c(n72-n73))
write.table(data.frame(Description=c("","Percentage of half-yearly in-house lobbyist observations excluded due to dyads with only one half-year with non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n72-n73)/n71,2))
write.table(data.frame(Description=c("","Total percentage of half-yearly in-house lobbyist observations excluded:")),row.names=F,col.names=F,quote=F)
print(round(c(n71-n73)/n71,2))

write.table(data.frame(Description=c("","Overall number of yearly observations before filtering:")),row.names=F,col.names=F,quote=F)
print(n67a1)
write.table(data.frame(Description=c("","Number of yearly observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(n67a1-n67a2)
write.table(data.frame(Description=c("","Percentage of yearly observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n67a1-n67a2)/n67a1,2))
write.table(data.frame(Description=c("","Number of yearly observations excluded due to relationships with only one year of non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(c(n67a2-n67a3))
write.table(data.frame(Description=c("","Percentage of yearly observations excluded due to relationships with only one year of non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n67a2-n67a3)/n67a1,2))
write.table(data.frame(Description=c("","Total percentage of yearly observations excluded:")),row.names=F,col.names=F,quote=F)
print(round(c(n67a1-n67a3)/n67a1,2))

write.table(data.frame(Description=c("","Number of yearly contract lobbyist observations before filtering:")),row.names=F,col.names=F,quote=F)
print(n6a1)
write.table(data.frame(Description=c("","Number of yearly contract lobbyist observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(n6a1-n6a2)
write.table(data.frame(Description=c("","Percentage of yearly contract lobbyist observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n6a1-n6a2)/n6a1,2))
write.table(data.frame(Description=c("","Number of yearly contract lobbyist observations excluded due to relationships with only one year of non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(c(n6a2-n6a3))
write.table(data.frame(Description=c("","Percentage of yearly contract lobbyist observations excluded due to relationships with only one year of non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n6a2-n6a3)/n6a1,2))
write.table(data.frame(Description=c("","Total percentage of yearly contract lobbyist observations excluded:")),row.names=F,col.names=F,quote=F)
print(round(c(n6a1-n6a3)/n6a1,2))

write.table(data.frame(Description=c("","Number of yearly in-house lobbyist observations before filtering:")),row.names=F,col.names=F,quote=F)
print(n7a1)
write.table(data.frame(Description=c("","Number of yearly in-house lobbyist observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(n7a1-n7a2)
write.table(data.frame(Description=c("","Percentage of yearly in-house lobbyist observations excluded due to zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n7a1-n7a2)/n7a1,2))
write.table(data.frame(Description=c("","Number of yearly in-house lobbyist observations excluded due to relationships with only one year of non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(c(n7a2-n7a3))
write.table(data.frame(Description=c("","Percentage of yearly in-house lobbyist observations excluded due to relationships with only one year of non-zero expenditures:")),row.names=F,col.names=F,quote=F)
print(round(c(n7a2-n7a3)/n7a1,2))
write.table(data.frame(Description=c("","Total percentage of yearly in-house lobbyist observations excluded:")),row.names=F,col.names=F,quote=F)
print(round(c(n7a1-n7a3)/n7a1,2))
